#####################################################################
# Appendix F
#####################################################################

rm(list=ls())

library(Hmisc)
library(exactRankTests)
library(sensitivitymw)
library(sensitivitymv)
library(stargazer)
library(msm)

##############################################
# Read data and prepare basic output
##############################################

load("~/Dropbox/Crime Brazil/06_replication/001_rematch_crime.Rdata")

table(d$t_ind)
table(d_match$t_ind)

##############################################################
# Sensitivity Analysis step 2
##############################################################  

treatment2 = d_rematch$outcome_crime_stronghand_wave4[d_rematch$t_ind==1]
control2 = d_rematch$outcome_crime_stronghand_wave4[d_rematch$t_ind==0]

sens.analysis.signedrank = function(diff, Gamma){
  rk = rank(abs(diff))
  s1 = 1*(diff>0)
  s2 = 1*(diff<0)
  W = sum(s1*rk)
  Eplus = sum((s1+s2)*rk*Gamma)/(1+Gamma)
  Eminus = sum((s1+s2)*rk)/(1+Gamma)
  V = sum((s1+s2)*rk*rk*Gamma)/((1+Gamma)^2)
  Dplus = (W-Eplus)/sqrt(V)
  Dminus = (W-Eminus)/sqrt(V)
  list(lowerbound = 1-pnorm(Dminus), upperbound = 1-pnorm(Dplus))
}

diff2 = treatment2-control2

gamma.100_s2 = round(sens.analysis.signedrank(diff2, 1.00)$upperbound, digits=4)
gamma.105_s2 = round(sens.analysis.signedrank(diff2, 1.05)$upperbound, digits=4) 
gamma.110_s2 = round(sens.analysis.signedrank(diff2, 1.10)$upperbound, digits=4) 
gamma.111_s2 = round(sens.analysis.signedrank(diff2, 1.11)$upperbound, digits=4) 
gamma.112_s2 = round(sens.analysis.signedrank(diff2, 1.12)$upperbound, digits=4)
gamma.113_s2 = round(sens.analysis.signedrank(diff2, 1.13)$upperbound, digits=4)
gamma.114_s2 = round(sens.analysis.signedrank(diff2, 1.14)$upperbound, digits=4)
gamma.115_s2 = round(sens.analysis.signedrank(diff2, 1.15)$upperbound, digits=4) 
gamma.116_s2 = round(sens.analysis.signedrank(diff2, 1.16)$upperbound, digits=4) 
gamma.117_s2 = round(sens.analysis.signedrank(diff2, 1.17)$upperbound, digits=4) 
gamma.118_s2 = round(sens.analysis.signedrank(diff2, 1.18)$upperbound, digits=4) 

gamma.100_s2
gamma.105_s2
gamma.110_s2
gamma.111_s2
gamma.112_s2
gamma.113_s2 
gamma.114_s2
gamma.115_s2
gamma.116_s2
gamma.117_s2 # non-sign
gamma.118_s2

# Amplification test
# amplify(gamma, lambda)
# lambda: odds of exposure to treatment (top)
# delta: odds of greater outcome (bottom)
amplify(1.16, c(2)) 
